stanEx4
options(scipen=100,digits=3)
library(cmdstanr)
options(mc.cores=parallel::detectCores())
options(cmdstanr_max_rows=1000)
library(gridExtra)execute mcmc sampling
goMCMC=function(mdl,data,smp=500,wrm=100,th=1){
mcmc=mdl$sample(
data=data,
seed=1,
chains=4,
iter_sampling=smp,
iter_warmup=wrm,
thin=th,
refresh=1000
)
mcmc
}see mcmc result and parameters
seeMCMC=function(mcmc,exc='',ch=T){ # exclude 'exc' parameters from seeing
print(mcmc)
prs=mcmc$metadata()$model_params[-1] # reject lp__
for(pr in prs){
if(grepl('^y',pr)) next # not show predictive value "y*" information
if(exc!='' && grepl(paste0('^',exc),pr)) next
drw=mcmc$draws(pr)
if(ch){
par(mfrow=c(4,1),mar=c(1,5,1,1))
drw[,1,] |> plot(type='l',xlab='',ylab=pr)
drw[,2,] |> plot(type='l',xlab='',ylab=pr)
drw[,3,] |> plot(type='l',xlab='',ylab=pr)
drw[,4,] |> plot(type='l',xlab='',ylab=pr)
par(mar=c(3,5,3,3))
}
par(mfrow=c(1,2))
drw |> hist(main=pr,xlab='')
drw |> density() |> plot(main=pr)
}
}fn=function(mdl,data,smp=500,wrm=100){
mle=mdl$optimize(data=data)
print(mle)
mcmc=goMCMC(mdl,data,smp,wrm)
mcmc$metadata()$stan_variables
seeMCMC(mcmc,'m')
y0=mcmc$draws('y0')
smy0=summary(y0)
grid.arrange(
qplot(y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))
par(mfrow=c(1,2))
hist(y-smy0$median,xlab='obs.-prd.',main='residual')
density(y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')
m1=mcmc$draws('m1')
smm1=summary(m1)
y1=mcmc$draws('y1')
smy1=summary(y1)
xy=tibble(x=x1,m=smm1$median,ml5=smm1$q5,mu5=smm1$q95,yl5=smy1$q5,yu5=smy1$q95)
qplot(x,y,col=I('red'))+
geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
geom_line(aes(x=x,y=m),xy,col='black')
}explanatory variable obs.x have noise
xN(x0.sx),yN(b0+b1*x0,s)
ex8-0.stan
normal regression without explanatory variable’s noise
data{
int N;
vector[N] x;
vector[N] y;
int N1;
vector[N1] x1;
}
parameters{
real b0;
real b1;
real<lower=0> s;
}
model{
y~normal(b0+b1*x,s);
}
generated quantities{
vector[N] m0;
vector[N] y0;
for(i in 1:N){
m0[i]=b0+b1*x[i];
y0[i]=normal_rng(m0[i],s);
}
vector[N1] m1;
vector[N1] y1;
for(i in 1:N1){
m1[i]=b0+b1*x1[i];
y1[i]=normal_rng(m1[i],s);
}
}ex9.stan
normal regression with explanatory variable’s noise
data{
int N;
vector[N] x;
vector[N] y;
int N1;
vector[N1] x10;
}
parameters{
real b0;
real b1;
real<lower=0> s;
real<lower=0> sx;
vector[N] x0;
}
model{
for(i in 1:N){
x[i]~normal(x0[i],sx);
y[i]~normal(b0+b1*x0[i],s);
}
}
generated quantities{
vector[N] m0;
vector[N] y0;
for(i in 1:N){
m0[i]=b0+b1*x0[i];
y0[i]=normal_rng(m0[i],s);
}
vector[N1] m1;
vector[N1] x1;
vector[N1] y1;
for(i in 1:N1){
x1[i]=normal_rng(x10[i],sx);
m1[i]=b0+b1*x10[i];
y1[i]=normal_rng(m1[i],s);
}
}n=20
x0=runif(n,0,20)
x=rnorm(n,x0,2)
y=rnorm(n,x0*2+5,2)
qplot(x,y)
n1=10
#explanatory variable do not has noise
x1=seq(min(x),max(x),length.out=n1) # new data fpr predcition
data=list(N=n,x=x,y=y,N1=n1,x1=x1)
mdl=cmdstan_model('./ex8-0.stan')
fn(mdl,data)## Initial log joint probability = -4384.8
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 24 -36.8162 0.000474379 0.00124184 0.9102 0.9102 43
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.2 seconds.
## variable estimate
## lp__ -36.82
## b0 6.43
## b1 1.75
## s 3.82
## m0[1] 28.17
## m0[2] 35.44
## m0[3] 10.55
## m0[4] 11.11
## m0[5] 14.44
## m0[6] 20.73
## m0[7] 43.22
## m0[8] 4.41
## m0[9] 35.08
## m0[10] 31.91
## m0[11] 37.97
## m0[12] 31.78
## m0[13] 29.86
## m0[14] 39.48
## m0[15] 30.89
## m0[16] 45.01
## m0[17] 16.60
## m0[18] 9.70
## m0[19] 4.38
## m0[20] 25.93
## y0[1] 23.02
## y0[2] 31.09
## y0[3] 14.82
## y0[4] 9.54
## y0[5] 15.89
## y0[6] 23.38
## y0[7] 41.94
## y0[8] 5.91
## y0[9] 33.80
## y0[10] 23.17
## y0[11] 42.86
## y0[12] 36.30
## y0[13] 32.92
## y0[14] 35.46
## y0[15] 36.00
## y0[16] 43.62
## y0[17] 23.45
## y0[18] 12.53
## y0[19] 7.67
## y0[20] 22.98
## m1[1] 4.38
## m1[2] 8.90
## m1[3] 13.41
## m1[4] 17.92
## m1[5] 22.44
## m1[6] 26.95
## m1[7] 31.47
## m1[8] 35.98
## m1[9] 40.49
## m1[10] 45.01
## y1[1] 0.92
## y1[2] 12.45
## y1[3] 9.97
## y1[4] 19.15
## y1[5] 25.82
## y1[6] 20.50
## y1[7] 32.22
## y1[8] 29.03
## y1[9] 43.35
## y1[10] 41.12
## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.3 seconds.
##
## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -36.99 -36.71 1.21 0.99 -39.39 -35.65 1.00 667 588
## b0 6.40 6.44 1.66 1.61 3.69 8.95 1.00 328 345
## b1 1.75 1.75 0.13 0.13 1.54 1.96 1.00 449 573
## s 4.31 4.24 0.76 0.75 3.26 5.59 1.01 1081 1294
## m0[1] 28.14 28.16 1.00 0.95 26.47 29.79 1.00 1544 1293
## m0[2] 35.40 35.42 1.25 1.21 33.33 37.41 1.00 2487 1410
## m0[3] 10.52 10.53 1.42 1.38 8.22 12.70 1.00 333 374
## m0[4] 11.08 11.10 1.39 1.36 8.84 13.24 1.00 335 391
## m0[5] 14.40 14.45 1.22 1.20 12.43 16.33 1.00 354 433
## m0[6] 20.70 20.71 1.00 0.99 19.07 22.35 1.00 494 578
## m0[7] 43.18 43.21 1.68 1.65 40.40 45.83 1.00 1376 1453
## m0[8] 4.38 4.43 1.78 1.77 1.46 7.13 1.00 328 348
## m0[9] 35.04 35.07 1.24 1.19 32.99 37.03 1.00 2524 1410
## m0[10] 31.87 31.89 1.10 1.08 30.04 33.68 1.00 2510 1526
## m0[11] 37.94 37.95 1.38 1.35 35.62 40.07 1.00 2055 1506
## m0[12] 31.75 31.76 1.10 1.07 29.93 33.55 1.00 2495 1526
## m0[13] 29.82 29.84 1.04 0.99 28.14 31.52 1.00 2069 1354
## m0[14] 39.44 39.46 1.46 1.43 36.99 41.72 1.00 1812 1369
## m0[15] 30.86 30.86 1.07 1.02 29.10 32.61 1.00 2346 1453
## m0[16] 44.97 45.01 1.80 1.75 42.01 47.79 1.00 1240 1342
## m0[17] 16.57 16.62 1.13 1.09 14.75 18.38 1.00 379 486
## m0[18] 9.67 9.69 1.46 1.42 7.28 11.93 1.00 331 364
## m0[19] 4.35 4.40 1.78 1.77 1.43 7.10 1.00 328 348
## m0[20] 25.89 25.90 0.96 0.93 24.31 27.51 1.00 1017 1183
## y0[1] 28.34 28.41 4.52 4.33 21.00 35.70 1.00 1819 1919
## y0[2] 35.33 35.46 4.70 4.55 27.18 43.14 1.00 2070 1765
## y0[3] 10.57 10.67 4.59 4.44 2.89 18.02 1.00 1506 1999
## y0[4] 10.95 10.95 4.56 4.25 3.35 18.43 1.00 1485 1643
## y0[5] 14.32 14.38 4.49 4.37 6.66 21.40 1.00 1445 1513
## y0[6] 20.78 20.83 4.38 4.26 13.42 27.86 1.00 2036 1972
## y0[7] 43.14 42.98 4.73 4.66 35.51 50.86 1.00 1825 1553
## y0[8] 4.49 4.51 4.83 4.71 -3.46 12.40 1.00 1247 1542
## y0[9] 34.89 34.81 4.49 4.33 27.56 42.38 1.00 2041 1804
## y0[10] 31.75 31.78 4.58 4.20 24.20 39.25 1.00 1966 2011
## y0[11] 37.97 37.98 4.51 4.33 30.69 45.34 1.00 2046 1856
## y0[12] 31.74 31.71 4.44 4.25 24.39 39.16 1.00 1857 1914
## y0[13] 29.79 29.82 4.45 4.32 22.35 37.12 1.00 2063 2056
## y0[14] 39.47 39.33 4.54 4.53 32.11 47.21 1.00 2055 1970
## y0[15] 30.78 30.81 4.48 4.38 23.36 38.00 1.00 2038 1819
## y0[16] 44.89 44.99 4.84 4.70 36.83 52.78 1.00 1866 1932
## y0[17] 16.55 16.49 4.55 4.34 8.99 23.82 1.00 1640 1918
## y0[18] 9.87 9.92 4.57 4.45 2.24 17.18 1.00 1341 1442
## y0[19] 4.45 4.46 4.67 4.63 -3.13 12.08 1.00 1133 1427
## y0[20] 25.98 25.95 4.55 4.35 18.50 33.57 1.00 2042 1927
## m1[1] 4.35 4.40 1.78 1.77 1.43 7.10 1.00 328 348
## m1[2] 8.87 8.90 1.51 1.47 6.40 11.21 1.00 330 369
## m1[3] 13.38 13.41 1.27 1.23 11.32 15.37 1.01 346 411
## m1[4] 17.89 17.92 1.08 1.06 16.15 19.62 1.00 404 537
## m1[5] 22.41 22.42 0.97 0.95 20.83 24.01 1.00 595 748
## m1[6] 26.92 26.93 0.98 0.94 25.31 28.55 1.00 1230 1196
## m1[7] 31.43 31.44 1.09 1.05 29.64 33.23 1.00 2447 1525
## m1[8] 35.94 35.97 1.28 1.24 33.82 37.98 1.00 2535 1393
## m1[9] 40.46 40.48 1.52 1.49 37.92 42.84 1.00 1672 1289
## m1[10] 44.97 45.01 1.80 1.75 42.01 47.79 1.00 1240 1342
## y1[1] 4.48 4.47 4.62 4.42 -2.85 12.06 1.00 1057 1380
## y1[2] 8.91 8.95 4.57 4.31 1.38 16.46 1.00 1584 1955
## y1[3] 13.43 13.34 4.59 4.37 6.12 20.94 1.00 1703 1876
## y1[4] 17.78 17.61 4.59 4.34 10.36 25.01 1.00 1565 1899
## y1[5] 22.57 22.61 4.54 4.33 15.06 29.91 1.00 1664 1855
## y1[6] 26.88 27.01 4.55 4.19 19.33 34.40 1.00 2017 2001
## y1[7] 31.36 31.41 4.55 4.27 24.03 38.85 1.00 2017 2083
## y1[8] 35.81 35.81 4.53 4.55 28.39 43.14 1.00 1943 1678
## y1[9] 40.44 40.38 4.66 4.52 32.78 48.02 1.00 1814 1917
## y1[10] 44.91 44.80 4.68 4.63 37.44 52.38 1.00 2113 1974
#explanatory variable has noise
x10=seq(min(x),max(x),length.out=n1) # new data fpr predcition
data=list(N=n,x=x,y=y,N1=n1,x10=x10)
mdl=cmdstan_model('./ex9.stan')
mcmc=goMCMC(mdl,data,wrm=500,smp=1000)## Running MCMC with 4 parallel chains...
##
## Chain 1 Iteration: 1 / 1500 [ 0%] (Warmup)
## Chain 1 Iteration: 501 / 1500 [ 33%] (Sampling)
## Chain 2 Iteration: 1 / 1500 [ 0%] (Warmup)
## Chain 3 Iteration: 1 / 1500 [ 0%] (Warmup)
## Chain 4 Iteration: 1 / 1500 [ 0%] (Warmup)
## Chain 2 Iteration: 501 / 1500 [ 33%] (Sampling)
## Chain 3 Iteration: 501 / 1500 [ 33%] (Sampling)
## Chain 4 Iteration: 501 / 1500 [ 33%] (Sampling)
## Chain 1 Iteration: 1500 / 1500 [100%] (Sampling)
## Chain 4 Iteration: 1500 / 1500 [100%] (Sampling)
## Chain 1 finished in 0.4 seconds.
## Chain 4 finished in 0.4 seconds.
## Chain 2 Iteration: 1500 / 1500 [100%] (Sampling)
## Chain 3 Iteration: 1500 / 1500 [100%] (Sampling)
## Chain 2 finished in 0.5 seconds.
## Chain 3 finished in 0.5 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.5 seconds.
## Total execution time: 0.6 seconds.seeMCMC(mcmc,'[m,x]')## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -41.21 -43.66 12.40 10.50 -57.54 -16.72 1.11 32 61
## b0 6.51 6.48 1.81 1.70 3.58 9.54 1.01 1065 1471
## b1 1.74 1.75 0.14 0.13 1.51 1.97 1.01 903 1072
## s 3.37 3.46 1.30 1.28 1.17 5.39 1.03 173 239
## sx 1.39 1.29 0.91 1.06 0.20 2.86 1.09 36 75
## x0[1] 13.98 13.49 1.76 1.67 11.90 17.09 1.03 179 1011
## x0[2] 16.54 16.55 1.06 0.75 14.77 18.28 1.02 2702 1516
## x0[3] 1.56 1.87 1.29 0.99 -0.88 3.19 1.01 458 1100
## x0[4] 3.81 3.46 1.46 1.35 1.98 6.46 1.02 174 1440
## x0[5] 4.52 4.54 1.06 0.76 2.77 6.24 1.04 2648 1775
## x0[6] 7.25 7.56 1.29 1.13 4.87 8.89 1.02 250 1150
## x0[7] 21.31 21.16 1.18 0.83 19.64 23.35 1.03 1465 948
## x0[8] -0.79 -0.93 1.16 0.86 -2.61 1.21 1.02 769 1578
## x0[9] 16.59 16.49 1.07 0.76 14.97 18.35 1.03 2167 1335
## x0[10] 13.70 13.95 1.26 1.14 11.41 15.35 1.02 364 1903
## x0[11] 18.80 18.51 1.31 1.04 17.17 21.13 1.01 466 973
## x0[12] 14.72 14.63 1.06 0.79 13.05 16.57 1.03 1784 1287
## x0[13] 13.70 13.60 1.07 0.83 12.05 15.52 1.02 1501 1435
## x0[14] 19.62 19.33 1.31 0.98 17.89 21.93 1.02 486 1111
## x0[15] 12.70 13.09 1.52 1.45 9.89 14.55 1.03 139 1158
## x0[16] 20.96 21.32 1.44 1.29 18.38 22.74 1.02 185 1288
## x0[17] 4.90 5.22 1.33 1.12 2.46 6.59 1.02 358 1201
## x0[18] 0.90 1.27 1.38 1.12 -1.70 2.55 1.02 365 983
## x0[19] -0.76 -0.92 1.14 0.89 -2.48 1.23 1.03 764 1527
## x0[20] 12.18 11.84 1.41 1.25 10.41 14.65 1.02 276 1413
## m0[1] 30.80 30.28 2.99 3.23 26.67 35.70 1.02 208 1579
## m0[2] 35.28 35.29 1.91 1.68 31.99 38.32 1.01 3888 2506
## m0[3] 9.28 9.38 2.27 2.32 5.62 12.90 1.02 343 1969
## m0[4] 13.16 12.89 2.64 2.82 9.22 17.58 1.02 235 1622
## m0[5] 14.40 14.40 1.91 1.67 11.30 17.51 1.01 3094 2649
## m0[6] 19.16 19.40 2.28 2.39 15.29 22.57 1.02 170 1706
## m0[7] 43.54 43.64 2.07 1.91 39.99 46.88 1.01 1784 2614
## m0[8] 5.19 5.23 2.20 2.08 1.57 8.73 1.00 1141 2147
## m0[9] 35.35 35.37 1.89 1.68 32.25 38.40 1.01 3345 2282
## m0[10] 30.34 30.51 2.27 2.41 26.55 33.70 1.02 276 2416
## m0[11] 39.19 39.10 2.26 2.31 35.67 42.86 1.01 512 2526
## m0[12] 32.12 32.13 1.88 1.60 28.97 35.16 1.01 2432 2278
## m0[13] 30.34 30.34 1.86 1.72 27.29 33.36 1.01 1756 2645
## m0[14] 40.60 40.58 2.22 2.26 36.96 44.15 1.01 518 2563
## m0[15] 28.61 29.11 2.71 2.93 24.01 32.39 1.03 105 1631
## m0[16] 42.95 43.09 2.66 2.99 38.65 47.03 1.02 182 2215
## m0[17] 15.08 15.27 2.34 2.53 11.23 18.63 1.02 310 2050
## m0[18] 8.13 8.26 2.39 2.56 4.27 12.00 1.02 338 1604
## m0[19] 5.24 5.26 2.18 2.07 1.65 8.65 1.01 865 1767
## m0[20] 27.69 27.36 2.43 2.55 24.11 31.77 1.02 303 2359
## y0[1] 30.75 31.15 4.70 4.76 22.72 37.48 1.01 425 2952
## y0[2] 35.33 35.30 4.11 3.69 28.53 42.22 1.00 4227 3618
## y0[3] 9.32 8.91 4.36 3.96 2.69 16.86 1.00 1659 2918
## y0[4] 13.10 13.50 4.48 4.34 5.27 19.71 1.01 693 2926
## y0[5] 14.41 14.34 4.16 3.53 7.71 21.36 1.00 3622 2945
## y0[6] 19.14 18.78 4.20 3.91 12.87 26.36 1.01 959 3038
## y0[7] 43.55 43.64 4.09 3.48 36.85 50.14 1.00 3382 2902
## y0[8] 5.14 5.36 4.20 3.75 -2.16 11.63 1.00 2523 2828
## y0[9] 35.29 35.38 4.03 3.53 28.42 41.86 1.00 3708 3299
## y0[10] 30.23 29.88 4.28 3.90 23.91 37.67 1.01 1501 3525
## y0[11] 39.18 39.49 4.18 3.76 31.81 45.52 1.00 1763 2431
## y0[12] 32.18 32.36 4.12 3.55 25.50 38.53 1.00 3823 3138
## y0[13] 30.38 30.47 4.01 3.60 23.49 36.64 1.00 4339 3596
## y0[14] 40.65 41.06 4.19 3.89 33.37 46.98 1.00 2191 3456
## y0[15] 28.60 28.10 4.52 4.47 22.19 36.61 1.01 664 2621
## y0[16] 42.88 42.35 4.53 4.39 36.17 50.94 1.01 886 3340
## y0[17] 15.05 14.72 4.25 3.93 8.64 22.68 1.00 1252 2241
## y0[18] 8.06 7.66 4.26 3.94 1.48 15.46 1.00 1402 3120
## y0[19] 5.28 5.49 4.18 3.70 -1.74 11.86 1.00 2682 2887
## y0[20] 27.62 28.11 4.25 4.09 20.13 33.75 1.01 884 3314
## m1[1] 4.47 4.42 1.95 1.84 1.32 7.70 1.01 1037 1453
## m1[2] 8.96 8.96 1.65 1.55 6.29 11.72 1.01 1108 1499
## m1[3] 13.45 13.45 1.38 1.30 11.24 15.76 1.01 1228 1564
## m1[4] 17.94 17.96 1.17 1.10 16.07 19.86 1.00 1437 1833
## m1[5] 22.43 22.46 1.05 1.00 20.75 24.13 1.00 1761 1962
## m1[6] 26.92 26.95 1.05 0.96 25.23 28.64 1.00 2009 2082
## m1[7] 31.42 31.44 1.17 1.08 29.54 33.23 1.00 1821 2022
## m1[8] 35.91 35.96 1.38 1.26 33.72 38.01 1.00 1464 1929
## m1[9] 40.40 40.46 1.64 1.51 37.74 42.92 1.00 1211 1494
## m1[10] 44.89 44.96 1.94 1.78 41.69 47.93 1.00 1096 1518
## x1[1] -1.18 -1.18 1.64 0.98 -3.91 1.55 1.04 4015 1415
## x1[2] 1.38 1.41 1.63 1.01 -1.43 4.15 1.04 4114 1970
## x1[3] 4.05 4.01 1.69 0.99 1.35 6.95 1.03 3549 1609
## x1[4] 6.61 6.59 1.63 1.00 3.93 9.36 1.03 3677 2227
## x1[5] 9.15 9.15 1.69 1.05 6.34 11.83 1.04 4111 2304
## x1[6] 11.70 11.70 1.64 0.96 8.99 14.38 1.03 3896 1756
## x1[7] 14.32 14.29 1.68 1.05 11.58 17.14 1.04 3890 2313
## x1[8] 16.92 16.91 1.63 1.00 14.20 19.61 1.04 4201 2090
## x1[9] 19.49 19.47 1.69 0.99 16.79 22.30 1.04 3808 2045
## x1[10] 22.02 22.05 1.67 0.99 19.24 24.76 1.04 3835 2267
## y1[1] 4.48 4.43 4.15 3.74 -2.13 11.30 1.00 2580 2745
## y1[2] 8.96 8.86 4.02 3.59 2.65 15.46 1.00 2525 3112
## y1[3] 13.49 13.47 3.84 3.34 7.36 19.88 1.00 3089 3219
## y1[4] 18.01 18.00 3.82 3.40 11.84 24.38 1.00 3381 3288
## y1[5] 22.49 22.42 3.76 3.13 16.32 28.69 1.01 3807 3476
## y1[6] 26.87 26.87 3.72 3.30 20.60 32.98 1.00 3802 2987
## y1[7] 31.39 31.40 3.79 3.27 25.27 37.70 1.00 3324 3011
## y1[8] 35.86 35.86 3.87 3.35 29.54 42.21 1.01 3019 2990
## y1[9] 40.39 40.42 3.97 3.48 33.89 46.91 1.01 3022 2701
## y1[10] 44.90 44.92 4.11 3.63 38.26 51.76 1.00 2757 2626
y0=mcmc$draws('y0')
smy0=summary(y0)
grid.arrange(
qplot(y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))
par(mfrow=c(1,2))
hist(y-smy0$median,xlab='obs.-prd.',main='residual')
density(y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')
m1=mcmc$draws('m1')
smm1=summary(m1)
x1=mcmc$draws('x1')
smx1=summary(x1)
y1=mcmc$draws('y1')
smy1=summary(y1)
xy=tibble(x=smx1$median,m=smm1$median,ml5=smm1$q5,mu5=smm1$q95,yl5=smy1$q5,yu5=smy1$q95)
qplot(x,y,col=I('red'))+
geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
geom_line(aes(x=x,y=m),xy,col='black')
outlier
ex10.stan
objective variable have outlier, y~cauchy(b0+b1*x,s)
data{
int N;
vector[N] x;
vector[N] y;
int N1;
vector[N1] x1;
}
parameters{
real b0;
real b1;
real<lower=0> s;
}
model{
y~cauchy(b0+b1*x,s);
}
generated quantities{
vector[N] m0;
vector[N] y0;
for(i in 1:N){
m0[i]=b0+b1*x[i];
y0[i]=cauchy_rng(m0[i],s);
}
vector[N1] m1;
vector[N1] y1;
for(i in 1:N1){
m1[i]=b0+b1*x1[i];
y1[i]=cauchy_rng(m1[i],s);
}
}n=20
x=runif(n,0,9)
y=rnorm(n,x*2+5,1)
x[1]=3
y[1]=25
qplot(x,y)
n1=10
x1=seq(min(x),max(x),length.out=n1) # new data fpr predcition
data=list(N=n,x=x,y=y,N1=n1,x1=x1)
mdl=cmdstan_model('./ex8-0.stan')
fn(mdl,data)## Initial log joint probability = -34076.2
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 24 -32.6737 0.000167123 0.000609913 0.8415 0.8415 58
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.
## variable estimate
## lp__ -32.67
## b0 6.72
## b1 1.82
## s 3.11
## m0[1] 12.17
## m0[2] 22.76
## m0[3] 10.33
## m0[4] 18.05
## m0[5] 16.65
## m0[6] 9.37
## m0[7] 19.72
## m0[8] 12.30
## m0[9] 11.54
## m0[10] 20.39
## m0[11] 8.59
## m0[12] 22.16
## m0[13] 17.16
## m0[14] 19.07
## m0[15] 8.69
## m0[16] 20.91
## m0[17] 14.00
## m0[18] 12.49
## m0[19] 15.69
## m0[20] 9.98
## y0[1] 15.32
## y0[2] 25.10
## y0[3] 9.12
## y0[4] 19.79
## y0[5] 18.27
## y0[6] 8.03
## y0[7] 16.79
## y0[8] 6.34
## y0[9] 13.10
## y0[10] 22.78
## y0[11] 6.79
## y0[12] 23.02
## y0[13] 13.50
## y0[14] 15.23
## y0[15] 12.08
## y0[16] 21.69
## y0[17] 17.13
## y0[18] 13.76
## y0[19] 12.63
## y0[20] 12.74
## m1[1] 8.59
## m1[2] 10.16
## m1[3] 11.74
## m1[4] 13.31
## m1[5] 14.89
## m1[6] 16.46
## m1[7] 18.04
## m1[8] 19.61
## m1[9] 21.19
## m1[10] 22.76
## y1[1] 10.11
## y1[2] 7.40
## y1[3] 11.55
## y1[4] 9.03
## y1[5] 16.79
## y1[6] 16.39
## y1[7] 15.51
## y1[8] 15.82
## y1[9] 24.80
## y1[10] 23.42
## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.3 seconds.
##
## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -33.14 -32.79 1.31 1.07 -35.64 -31.74 1.01 493 637
## b0 6.66 6.74 1.60 1.57 3.94 9.11 1.01 311 301
## b1 1.83 1.82 0.30 0.29 1.36 2.32 1.01 406 472
## s 3.52 3.42 0.66 0.63 2.65 4.75 1.00 1236 1347
## m0[1] 12.15 12.18 0.95 0.92 10.60 13.64 1.01 389 423
## m0[2] 22.80 22.81 1.49 1.44 20.34 25.28 1.00 1075 1178
## m0[3] 10.30 10.35 1.13 1.11 8.38 12.08 1.01 331 346
## m0[4] 18.07 18.07 0.95 0.89 16.44 19.60 1.00 2043 1281
## m0[5] 16.65 16.65 0.85 0.80 15.22 18.01 1.00 1644 1248
## m0[6] 9.33 9.38 1.25 1.21 7.23 11.27 1.01 319 322
## m0[7] 19.74 19.74 1.11 1.04 17.88 21.56 1.00 1655 1353
## m0[8] 12.27 12.29 0.94 0.91 10.75 13.78 1.01 395 454
## m0[9] 11.51 11.55 1.01 0.98 9.85 13.11 1.01 360 381
## m0[10] 20.42 20.42 1.19 1.14 18.45 22.35 1.00 1462 1247
## m0[11] 8.54 8.61 1.35 1.32 6.25 10.60 1.01 314 305
## m0[12] 22.20 22.20 1.41 1.37 19.89 24.52 1.00 1148 1121
## m0[13] 17.17 17.17 0.88 0.82 15.67 18.57 1.00 1884 1298
## m0[14] 19.09 19.08 1.04 0.98 17.34 20.79 1.00 1850 1281
## m0[15] 8.64 8.71 1.33 1.30 6.37 10.68 1.01 314 310
## m0[16] 20.94 20.94 1.25 1.20 18.87 22.97 1.00 1351 1240
## m0[17] 13.99 13.99 0.84 0.80 12.63 15.33 1.00 574 886
## m0[18] 12.47 12.48 0.92 0.91 10.96 13.95 1.01 407 449
## m0[19] 15.69 15.69 0.82 0.77 14.30 17.00 1.00 1152 1186
## m0[20] 9.94 10.00 1.17 1.15 7.95 11.77 1.01 326 324
## y0[1] 12.14 12.18 3.77 3.51 5.96 18.32 1.00 1650 1893
## y0[2] 22.72 22.82 3.83 3.62 16.35 29.05 1.00 1907 1881
## y0[3] 10.23 10.21 3.76 3.47 4.11 16.37 1.00 1314 1609
## y0[4] 18.21 18.20 3.67 3.57 12.41 24.24 1.00 1961 1703
## y0[5] 16.64 16.61 3.76 3.77 10.73 22.71 1.00 2086 1822
## y0[6] 9.49 9.51 3.79 3.67 3.22 15.53 1.00 1450 1366
## y0[7] 19.64 19.63 3.78 3.50 13.51 25.80 1.00 2064 1966
## y0[8] 12.12 12.10 3.78 3.54 6.11 18.39 1.00 1785 1905
## y0[9] 11.55 11.47 3.74 3.56 5.66 17.51 1.00 1579 1643
## y0[10] 20.40 20.38 3.86 3.64 14.17 26.61 1.00 1740 1849
## y0[11] 8.55 8.54 3.85 3.71 2.14 14.77 1.00 1139 1591
## y0[12] 22.12 22.08 3.77 3.61 15.94 28.49 1.00 1727 1961
## y0[13] 17.17 17.20 3.81 3.60 10.52 23.18 1.00 1988 1971
## y0[14] 19.14 19.13 3.77 3.69 13.03 25.42 1.00 2137 1822
## y0[15] 8.68 8.73 3.97 3.94 2.02 15.29 1.00 1382 1380
## y0[16] 20.95 20.90 3.90 3.79 14.60 27.54 1.00 2005 1633
## y0[17] 14.00 13.85 3.73 3.38 8.00 20.19 1.00 1844 1807
## y0[18] 12.54 12.44 3.83 3.59 6.35 18.76 1.00 1689 1966
## y0[19] 15.60 15.66 3.62 3.39 9.70 21.49 1.00 1991 1658
## y0[20] 9.93 10.02 3.71 3.45 3.69 15.84 1.00 1340 1697
## m1[1] 8.54 8.61 1.35 1.32 6.25 10.60 1.01 314 305
## m1[2] 10.13 10.18 1.15 1.13 8.18 11.93 1.01 329 336
## m1[3] 11.71 11.74 0.99 0.97 10.09 13.28 1.01 368 390
## m1[4] 13.30 13.30 0.87 0.85 11.87 14.68 1.01 478 613
## m1[5] 14.88 14.88 0.82 0.78 13.53 16.20 1.00 789 1107
## m1[6] 16.46 16.46 0.85 0.79 15.04 17.80 1.00 1549 1298
## m1[7] 18.05 18.05 0.95 0.89 16.43 19.58 1.00 2044 1281
## m1[8] 19.63 19.63 1.10 1.03 17.79 21.44 1.00 1687 1353
## m1[9] 21.22 21.21 1.29 1.24 19.11 23.30 1.00 1298 1239
## m1[10] 22.80 22.81 1.49 1.44 20.34 25.28 1.00 1075 1178
## y1[1] 8.73 8.82 3.78 3.58 2.54 14.87 1.00 1159 1515
## y1[2] 10.15 10.17 3.76 3.57 3.95 16.15 1.00 1435 1580
## y1[3] 11.64 11.59 3.75 3.65 5.63 17.58 1.00 1866 1842
## y1[4] 13.24 13.19 3.64 3.43 7.33 19.35 1.00 1619 1582
## y1[5] 14.94 14.93 3.69 3.63 8.94 20.96 1.00 1832 1763
## y1[6] 16.46 16.35 3.68 3.51 10.59 22.41 1.00 1912 1854
## y1[7] 17.97 17.95 3.61 3.37 11.94 23.97 1.00 1896 1806
## y1[8] 19.66 19.67 3.73 3.58 13.51 25.84 1.00 1849 1861
## y1[9] 21.12 21.06 3.91 3.56 14.67 27.45 1.00 1744 1675
## y1[10] 22.81 22.82 3.88 3.58 16.53 29.10 1.00 1383 1582
mdl=cmdstan_model('./ex10.stan')
fn(mdl,data)## Initial log joint probability = -100.194
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 21 -13.2521 1.9685e-05 0.00237116 0.3438 0.3438 28
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.
## variable estimate
## lp__ -13.25
## b0 5.49
## b1 1.95
## s 0.58
## m0[1] 11.33
## m0[2] 22.65
## m0[3] 9.36
## m0[4] 17.62
## m0[5] 16.11
## m0[6] 8.32
## m0[7] 19.40
## m0[8] 11.46
## m0[9] 10.65
## m0[10] 20.12
## m0[11] 7.49
## m0[12] 22.01
## m0[13] 16.66
## m0[14] 18.70
## m0[15] 7.60
## m0[16] 20.67
## m0[17] 13.28
## m0[18] 11.67
## m0[19] 15.09
## m0[20] 8.98
## y0[1] 10.98
## y0[2] 23.16
## y0[3] 8.17
## y0[4] 10.04
## y0[5] 17.18
## y0[6] 7.49
## y0[7] 20.39
## y0[8] 11.94
## y0[9] 11.27
## y0[10] 20.19
## y0[11] 0.69
## y0[12] 22.09
## y0[13] 16.35
## y0[14] 25.24
## y0[15] 8.13
## y0[16] 21.04
## y0[17] 13.56
## y0[18] 11.65
## y0[19] 13.20
## y0[20] 9.77
## m1[1] 7.49
## m1[2] 9.17
## m1[3] 10.86
## m1[4] 12.54
## m1[5] 14.23
## m1[6] 15.91
## m1[7] 17.60
## m1[8] 19.28
## m1[9] 20.97
## m1[10] 22.65
## y1[1] 7.18
## y1[2] 9.58
## y1[3] 3.83
## y1[4] 12.06
## y1[5] 13.83
## y1[6] 16.78
## y1[7] 17.07
## y1[8] 19.85
## y1[9] 20.83
## y1[10] 22.43
## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.
##
## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -15.55 -15.21 1.39 1.16 -18.21 -13.99 1.00 625 862
## b0 5.46 5.46 0.51 0.49 4.63 6.29 1.01 494 555
## b1 1.95 1.95 0.09 0.09 1.80 2.09 1.00 592 884
## s 0.76 0.73 0.25 0.24 0.40 1.21 1.00 598 386
## m0[1] 11.30 11.31 0.30 0.29 10.79 11.79 1.01 586 675
## m0[2] 22.63 22.63 0.43 0.37 21.91 23.32 1.00 1168 1134
## m0[3] 9.33 9.34 0.36 0.33 8.74 9.92 1.01 522 622
## m0[4] 17.59 17.60 0.28 0.25 17.13 18.03 1.00 1834 1242
## m0[5] 16.09 16.09 0.25 0.23 15.66 16.48 1.00 1483 1316
## m0[6] 8.30 8.30 0.40 0.37 7.64 8.93 1.01 506 574
## m0[7] 19.38 19.38 0.32 0.28 18.84 19.87 1.00 1630 1346
## m0[8] 11.43 11.44 0.30 0.28 10.93 11.91 1.01 594 675
## m0[9] 10.62 10.63 0.32 0.30 10.09 11.14 1.01 556 634
## m0[10] 20.10 20.10 0.34 0.30 19.53 20.62 1.00 1512 1265
## m0[11] 7.46 7.47 0.43 0.40 6.75 8.16 1.01 499 565
## m0[12] 21.99 21.99 0.41 0.35 21.31 22.63 1.00 1238 1107
## m0[13] 16.64 16.64 0.26 0.23 16.19 17.05 1.00 1647 1299
## m0[14] 18.68 18.68 0.30 0.26 18.18 19.14 1.00 1768 1333
## m0[15] 7.57 7.57 0.42 0.39 6.86 8.25 1.01 499 561
## m0[16] 20.65 20.65 0.36 0.31 20.05 21.21 1.00 1432 1195
## m0[17] 13.25 13.26 0.26 0.24 12.81 13.68 1.01 767 811
## m0[18] 11.64 11.65 0.29 0.28 11.15 12.12 1.01 608 675
## m0[19] 15.06 15.06 0.25 0.23 14.63 15.46 1.00 1178 1217
## m0[20] 8.95 8.96 0.37 0.35 8.33 9.56 1.01 516 635
## y0[1] 11.29 11.29 15.01 1.20 5.85 15.68 1.00 2057 1942
## y0[2] 22.45 22.62 21.64 1.14 17.93 27.19 1.00 1936 1883
## y0[3] 10.11 9.38 75.82 1.19 4.71 14.77 1.00 1873 1965
## y0[4] 17.72 17.55 17.06 1.21 11.91 22.59 1.00 1878 1802
## y0[5] 17.10 16.11 25.82 1.13 11.40 21.05 1.00 1809 1753
## y0[6] 8.57 8.29 16.23 1.16 3.45 12.97 1.00 1933 1655
## y0[7] 20.95 19.35 47.55 1.15 14.78 24.41 1.00 1748 1545
## y0[8] 7.32 11.45 157.92 1.19 6.45 16.60 1.00 1915 1894
## y0[9] 7.97 10.61 303.88 1.22 5.28 15.64 1.00 1890 1799
## y0[10] 19.99 20.03 23.97 1.17 14.84 24.56 1.00 1715 1774
## y0[11] 7.41 7.51 12.27 1.22 3.15 12.09 1.00 1769 2097
## y0[12] 22.19 21.96 20.59 1.15 17.69 26.58 1.00 1991 1929
## y0[13] 24.35 16.64 326.61 1.13 12.65 21.32 1.00 2140 1933
## y0[14] 18.95 18.68 27.92 1.14 13.77 23.91 1.00 2059 1894
## y0[15] 7.90 7.54 79.37 1.24 1.85 12.27 1.00 1636 1905
## y0[16] 21.45 20.66 17.26 1.13 16.37 24.86 1.00 1874 1823
## y0[17] 13.66 13.22 35.95 1.19 7.68 17.98 1.00 1807 1643
## y0[18] 11.21 11.59 20.90 1.11 6.38 15.87 1.00 2180 2037
## y0[19] 15.40 15.06 22.02 1.14 10.38 20.18 1.00 1978 1865
## y0[20] 7.02 8.95 78.36 1.18 3.45 13.90 1.00 1907 1894
## m1[1] 7.46 7.47 0.43 0.40 6.75 8.16 1.01 499 565
## m1[2] 9.15 9.15 0.37 0.34 8.55 9.74 1.01 519 635
## m1[3] 10.83 10.84 0.31 0.30 10.31 11.34 1.01 564 647
## m1[4] 12.52 12.52 0.27 0.26 12.06 12.96 1.01 681 745
## m1[5] 14.20 14.21 0.25 0.24 13.78 14.62 1.00 942 954
## m1[6] 15.89 15.89 0.25 0.23 15.47 16.28 1.00 1433 1256
## m1[7] 17.57 17.58 0.27 0.25 17.11 18.01 1.00 1833 1242
## m1[8] 19.26 19.26 0.32 0.27 18.73 19.75 1.00 1653 1346
## m1[9] 20.94 20.95 0.37 0.31 20.33 21.53 1.00 1388 1135
## m1[10] 22.63 22.63 0.43 0.37 21.91 23.32 1.00 1168 1134
## y1[1] 5.10 7.47 163.26 1.24 3.06 12.61 1.01 1460 1997
## y1[2] 7.92 9.13 51.20 1.14 5.34 13.25 1.00 1651 1835
## y1[3] 10.32 10.84 12.90 1.09 6.68 15.52 1.00 1999 1908
## y1[4] 12.98 12.53 25.21 1.13 7.52 17.34 1.00 1974 1846
## y1[5] 13.85 14.20 10.99 1.11 9.00 18.71 1.00 2036 1959
## y1[6] 16.96 15.89 42.38 1.16 12.01 20.33 1.00 1904 1974
## y1[7] 17.51 17.58 83.77 1.17 12.98 21.73 1.00 1925 1792
## y1[8] 19.24 19.27 16.73 1.08 14.91 23.52 1.00 1922 1671
## y1[9] 14.26 20.95 328.50 1.20 16.27 25.62 1.00 1944 1845
## y1[10] 23.20 22.64 24.25 1.25 17.89 28.18 1.00 2077 1872
censored data
objective variable has NA
(x,y) i=1-N
(x0,y0) i=1-N0
x1 i=1-N1, y1=NA
(x,y)~N((mx,my),(sx2,sy2,sxy))
(x0,y0)~N((mx,my),(sx2,sy2,sxy))
x1~N(mx,sx2)ex11-0.stan
data{
int N0;
array[N0] vector[2] xy;
int N1;
vector[N1] x1;
}
parameters{
vector[2] m;
cov_matrix[2] s;
}
model{
target+=multi_normal_lpdf(xy | m, s);
x1~normal(m[1],s[1,1]^.5);
}
generated quantities{
vector[2] xy1;
xy1=multi_normal_rng(m,s);
real cr;
cr=s[1,2]/(s[1,1]*s[2,2])^.5;
}n=30
x=runif(n,0,9)
y=rnorm(n,10+3*x,4)
cor(x,y)## [1] 0.889qplot(x,y)
L=4
n0=sum(x>L)
x0=x[x>L]
y0=y[x>L]
x1=x[x<=L]
n1=sum(x<=L)
cor(x0,y0)## [1] 0.789qplot(x0,y0)
mdl=cmdstan_model('./ex11-0.stan')
data=list(N0=n0,xy=matrix(c(x0,y0),ncol=2),N1=n1,x1=x1)
mle=mdl$optimize(data=data)## Initial log joint probability = -1794.89
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 49 -127.4 0.000806172 0.00104209 0.8332 0.8332 70
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -127.40
## m[1] 5.73
## m[2] 27.80
## s[1,1] 5.00
## s[2,1] 17.09
## s[1,2] 17.09
## s[2,2] 73.08
## xy1[1] 6.13
## xy1[2] 27.00
## cr 0.89mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.3 seconds.
## Chain 2 finished in 0.4 seconds.
## Chain 3 finished in 0.3 seconds.
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 finished in 0.4 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.3 seconds.
## Total execution time: 0.5 seconds.seeMCMC(mcmc,ch=F)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -123.59 -123.24 1.75 1.52 -126.99 -121.47 1.01 543 652
## m[1] 5.74 5.74 0.46 0.47 5.00 6.50 1.01 573 639
## m[2] 27.84 27.87 1.89 1.93 24.75 30.85 1.01 506 599
## s[1,1] 6.63 6.24 2.11 1.74 4.06 10.40 1.01 858 840
## s[2,1] 22.89 21.49 8.70 7.19 12.06 38.83 1.01 526 516
## s[1,2] 22.89 21.49 8.70 7.19 12.06 38.83 1.01 526 516
## s[2,2] 101.75 94.20 41.53 33.51 51.12 179.53 1.02 523 523
## xy1[1] 5.82 5.76 2.51 2.47 1.65 9.96 1.00 1864 1745
## xy1[2] 28.02 27.99 9.82 9.31 11.46 44.56 1.00 1858 1758
## cr 0.88 0.89 0.07 0.05 0.76 0.95 1.00 617 866
xy=mcmc$draws('xy1')
cor(xy[,,1],xy[,,2])## [1] 0.871qplot(xy[,,1],xy[,,2])
objective variable is censored
y i=1-N, y~N(m,s)
actual ya i=1-Na
lower censored yl i=1-Nl, y<L, P(y<L)=cdf(L-m /s)
upper censored yu i=1-Nu, y>U, P(y>U)=ccdf(U-m /s)
cdf(z) cumulative normal density function, P((-Infinity,z],z~N(0,1))
ccdf(z) complementary CDF, P([z,Infinity),z~N(0,1))
P(y | x,m,s)=P(ya i=1-Na)* P(yl i=1-Nl)* P(yu i=1-Nu)ex11-1.stan
data{
int N;
vector[N] x;
vector[N] y;
real L;
int Nl;
vector[Nl] xl;
real U;
int Nu;
vector[Nu] xu;
int N1;
vector[N1] x1;
}
parameters{
real b0;
real b1;
real<lower=0> s;
}
model{
y~normal(b0+b1*x,s);
for(i in 1:Nl)
target+=normal_lcdf(L | b0+b1*xl[i],s);
for(i in 1:Nu)
target+=normal_lccdf(U | b0+b1*xu[i],s);
}
generated quantities{
vector[N] m0;
vector[N] y0;
for(i in 1:N){
m0[i]=b0+b1*x[i];
y0[i]=normal_rng(m0[i],s);
}
vector[N1] m1;
vector[N1] y1;
for(i in 1:N1){
m1[i]=b0+b1*x1[i];
y1[i]=normal_rng(m1[i],s);
}
}n0=20
x=runif(n0,0,9)
y=rnorm(n0,10+3*x,4)
L=15
y[y<L]=L
nl=sum(y==L)
U=30
y[y>U]=U
nu=sum(y==U)
n=n0-nl-nu
qplot(x,y)
xy0=tibble(x=x,y=y)
xya=filter(xy0, y>L & y<U)
xyl=filter(xy0, y==L)
xyu=filter(xy0, y==U)
n1=10
x1=seq(min(x),max(x),length.out=n1) # new data fpr predcitionmdl=cmdstan_model('./ex8-0.stan')
data=list(N=n,x=xya$x,y=xya$y,N1=n1,x1=x1)
mle=mdl$optimize(data=data)## Initial log joint probability = -566.386
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 29 -10.3006 0.000345667 0.000330398 0.9893 0.9893 41
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -10.30
## b0 13.44
## b1 2.02
## s 1.55
## m0[1] 25.98
## m0[2] 23.61
## m0[3] 21.55
## m0[4] 19.65
## m0[5] 16.19
## m0[6] 21.79
## m0[7] 23.68
## m0[8] 16.39
## m0[9] 28.83
## m0[10] 26.16
## m0[11] 26.50
## y0[1] 23.60
## y0[2] 21.30
## y0[3] 22.02
## y0[4] 22.74
## y0[5] 13.60
## y0[6] 21.89
## y0[7] 23.99
## y0[8] 17.77
## y0[9] 29.19
## y0[10] 24.87
## y0[11] 25.03
## m1[1] 13.55
## m1[2] 15.47
## m1[3] 17.40
## m1[4] 19.33
## m1[5] 21.26
## m1[6] 23.18
## m1[7] 25.11
## m1[8] 27.04
## m1[9] 28.97
## m1[10] 30.89
## y1[1] 12.20
## y1[2] 16.28
## y1[3] 18.75
## y1[4] 18.52
## y1[5] 19.20
## y1[6] 21.97
## y1[7] 26.35
## y1[8] 27.96
## y1[9] 29.76
## y1[10] 33.49mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.seeMCMC(mcmc,'m')## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -11.70 -11.29 1.56 1.23 -14.76 -10.05 1.01 529 632
## b0 13.27 13.30 1.70 1.57 10.45 15.86 1.02 319 318
## b1 2.06 2.05 0.34 0.32 1.53 2.62 1.01 361 388
## s 2.03 1.90 0.64 0.48 1.30 3.29 1.00 726 905
## m0[1] 26.03 26.00 0.86 0.72 24.68 27.37 1.00 1448 1373
## m0[2] 23.62 23.59 0.68 0.58 22.58 24.70 1.00 1683 1288
## m0[3] 21.52 21.52 0.70 0.64 20.44 22.61 1.00 663 846
## m0[4] 19.59 19.60 0.84 0.76 18.22 20.92 1.01 405 491
## m0[5] 16.06 16.11 1.29 1.18 13.97 18.06 1.02 326 328
## m0[6] 21.76 21.76 0.69 0.62 20.69 22.84 1.00 738 872
## m0[7] 23.68 23.65 0.68 0.58 22.64 24.77 1.00 1712 1272
## m0[8] 16.27 16.31 1.26 1.16 14.21 18.22 1.02 327 328
## m0[9] 28.92 28.90 1.22 1.02 26.99 30.95 1.00 733 894
## m0[10] 26.21 26.18 0.88 0.74 24.85 27.59 1.00 1369 1373
## m0[11] 26.55 26.52 0.92 0.78 25.14 27.99 1.00 1258 1246
## y0[1] 26.01 26.02 2.33 2.11 22.14 29.66 1.00 1879 1933
## y0[2] 23.62 23.66 2.28 2.03 19.96 27.17 1.00 1699 1843
## y0[3] 21.54 21.55 2.25 2.07 17.90 25.16 1.00 1402 1641
## y0[4] 19.45 19.52 2.34 2.11 15.78 23.05 1.00 1007 1165
## y0[5] 16.04 16.13 2.49 2.30 11.95 19.97 1.00 765 925
## y0[6] 21.66 21.69 2.28 1.97 17.95 25.15 1.00 1560 1595
## y0[7] 23.66 23.67 2.25 1.98 20.12 27.27 1.00 1847 1916
## y0[8] 16.31 16.38 2.46 2.22 12.23 20.31 1.00 769 1233
## y0[9] 28.88 28.83 2.40 2.16 25.14 32.80 1.00 1453 1571
## y0[10] 26.20 26.18 2.28 2.01 22.49 29.89 1.00 1745 1767
## y0[11] 26.58 26.47 2.34 2.10 22.95 30.40 1.00 1695 1778
## m1[1] 13.38 13.41 1.68 1.55 10.59 15.95 1.02 319 318
## m1[2] 15.34 15.39 1.39 1.29 13.06 17.48 1.02 323 322
## m1[3] 17.30 17.32 1.12 1.02 15.46 19.04 1.02 336 338
## m1[4] 19.26 19.28 0.88 0.80 17.82 20.65 1.01 387 468
## m1[5] 21.22 21.23 0.71 0.64 20.11 22.33 1.00 592 702
## m1[6] 23.18 23.15 0.67 0.57 22.16 24.24 1.00 1455 1271
## m1[7] 25.14 25.11 0.78 0.66 23.93 26.35 1.00 1755 1379
## m1[8] 27.10 27.07 0.98 0.83 25.58 28.64 1.00 1060 1143
## m1[9] 29.06 29.04 1.24 1.04 27.09 31.14 1.00 719 885
## m1[10] 31.02 31.01 1.52 1.33 28.60 33.51 1.00 591 740
## y1[1] 13.42 13.44 2.74 2.44 8.82 17.67 1.01 562 766
## y1[2] 15.31 15.37 2.61 2.37 11.08 19.36 1.00 862 1079
## y1[3] 17.29 17.30 2.42 2.14 13.43 21.06 1.00 992 1316
## y1[4] 19.26 19.33 2.31 2.09 15.50 22.88 1.00 1349 1469
## y1[5] 21.21 21.21 2.29 2.04 17.58 25.02 1.00 1822 1652
## y1[6] 23.17 23.11 2.24 2.02 19.58 26.73 1.00 1748 1758
## y1[7] 25.11 25.00 2.23 2.02 21.61 28.61 1.00 1903 1688
## y1[8] 27.09 27.07 2.30 2.08 23.44 30.77 1.00 1687 1605
## y1[9] 29.07 28.94 2.48 2.28 25.17 33.19 1.00 1304 1257
## y1[10] 31.01 30.95 2.62 2.44 26.94 35.23 1.00 997 1316
y0=mcmc$draws('y0')
smy0=summary(y0)
grid.arrange(
qplot(xya$y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))
par(mfrow=c(1,2))
hist(xya$y-smy0$median,xlab='obs.-prd.',main='residual')
density(xya$y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')
m1=mcmc$draws('m1')
smm1=summary(m1)
y1=mcmc$draws('y1')
smy1=summary(y1)
xy=tibble(x=x1,m=smm1$median,ml5=smm1$q5,mu5=smm1$q95,yl5=smy1$q5,yu5=smy1$q95)
qplot(x,y,col=I('red'))+
geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
geom_line(aes(x=x,y=m),xy,col='black')
data=list(N=n,x=xya$x,y=xya$y,
L=L,Nl=nl,xl=xyl$x,
U=U,Nu=nu,xu=xyu$x,
N1=n1,x1=x1)
mdl=cmdstan_model('./ex11-1.stan')
mle=mdl$optimize(data=data)## Rejecting initial value:
## Log probability evaluates to log(0), i.e. negative infinity.
## Stan can't start sampling from this initial value.
## Initial log joint probability = -52.8227
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 23 -23.4739 0.00204867 0.000134233 1 1 26
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -23.47
## b0 8.42
## b1 2.97
## s 3.14
## m0[1] 26.84
## m0[2] 23.36
## m0[3] 20.33
## m0[4] 17.54
## m0[5] 12.45
## m0[6] 20.69
## m0[7] 23.46
## m0[8] 12.75
## m0[9] 31.02
## m0[10] 27.10
## m0[11] 27.60
## y0[1] 27.44
## y0[2] 20.73
## y0[3] 18.51
## y0[4] 12.20
## y0[5] 16.59
## y0[6] 24.23
## y0[7] 24.74
## y0[8] 17.58
## y0[9] 31.11
## y0[10] 29.86
## y0[11] 25.67
## m1[1] 8.57
## m1[2] 11.40
## m1[3] 14.23
## m1[4] 17.07
## m1[5] 19.90
## m1[6] 22.73
## m1[7] 25.56
## m1[8] 28.39
## m1[9] 31.23
## m1[10] 34.06
## y1[1] 9.63
## y1[2] 15.14
## y1[3] 16.92
## y1[4] 15.03
## y1[5] 22.33
## y1[6] 27.32
## y1[7] 24.73
## y1[8] 30.02
## y1[9] 33.52
## y1[10] 32.69mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.seeMCMC(mcmc,'m',ch=T)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -24.01 -23.62 1.39 1.07 -26.58 -22.52 1.01 581 675
## b0 7.14 7.37 2.89 2.45 2.05 11.53 1.02 234 366
## b1 3.23 3.19 0.58 0.51 2.36 4.24 1.01 288 357
## s 4.16 3.93 1.27 0.98 2.66 6.56 1.01 718 751
## m0[1] 27.13 27.04 1.44 1.35 25.02 29.53 1.00 1047 804
## m0[2] 23.35 23.33 1.13 1.06 21.64 25.21 1.00 1387 1059
## m0[3] 20.07 20.14 1.15 1.04 18.15 21.83 1.01 534 724
## m0[4] 17.04 17.17 1.41 1.22 14.63 19.18 1.02 280 474
## m0[5] 11.51 11.71 2.18 1.84 7.77 14.86 1.02 222 360
## m0[6] 20.45 20.50 1.13 1.02 18.57 22.20 1.01 592 730
## m0[7] 23.46 23.43 1.13 1.07 21.75 25.33 1.00 1409 1096
## m0[8] 11.84 12.03 2.13 1.80 8.20 15.12 1.02 221 360
## m0[9] 31.67 31.52 2.06 1.88 28.65 35.11 1.00 617 618
## m0[10] 27.42 27.31 1.47 1.39 25.26 29.87 1.00 999 800
## m0[11] 27.96 27.85 1.54 1.43 25.71 30.54 1.00 917 800
## y0[1] 27.25 27.15 4.53 4.14 20.03 34.58 1.00 1953 1723
## y0[2] 23.15 23.22 4.47 3.97 15.96 30.09 1.00 2060 1726
## y0[3] 20.23 20.18 4.56 3.94 13.01 27.65 1.00 1891 1605
## y0[4] 17.06 17.26 4.58 3.93 9.15 24.12 1.00 1465 1417
## y0[5] 11.44 11.68 4.82 4.49 3.53 18.73 1.01 910 1239
## y0[6] 20.49 20.62 4.52 3.96 12.96 27.51 1.00 2020 1775
## y0[7] 23.47 23.54 4.49 4.08 16.03 30.55 1.00 1931 1457
## y0[8] 11.93 12.15 4.86 4.26 3.50 19.56 1.00 1053 1224
## y0[9] 31.80 31.56 4.87 4.40 24.22 39.83 1.00 1417 1480
## y0[10] 27.30 27.15 4.68 4.25 19.92 35.02 1.00 1884 1950
## y0[11] 27.88 27.83 4.77 4.33 20.48 35.38 1.00 1939 1635
## m1[1] 7.30 7.54 2.86 2.43 2.26 11.66 1.02 233 365
## m1[2] 10.38 10.60 2.36 1.99 6.26 14.00 1.02 224 354
## m1[3] 13.45 13.61 1.89 1.61 10.17 16.32 1.02 222 389
## m1[4] 16.52 16.65 1.47 1.28 13.99 18.74 1.02 263 449
## m1[5] 19.60 19.68 1.18 1.07 17.63 21.40 1.01 474 662
## m1[6] 22.67 22.66 1.11 1.01 20.93 24.50 1.00 1200 1002
## m1[7] 25.74 25.68 1.29 1.24 23.77 27.90 1.00 1290 1070
## m1[8] 28.82 28.70 1.65 1.53 26.42 31.59 1.00 811 736
## m1[9] 31.89 31.75 2.10 1.89 28.86 35.39 1.00 606 624
## m1[10] 34.96 34.75 2.59 2.29 31.20 39.27 1.01 504 502
## y1[1] 7.27 7.51 5.21 4.73 -1.41 14.91 1.00 792 994
## y1[2] 10.27 10.45 5.02 4.41 2.07 17.79 1.01 945 1168
## y1[3] 13.32 13.53 4.71 4.29 5.51 20.61 1.00 1032 1147
## y1[4] 16.47 16.62 4.65 4.11 8.66 23.56 1.00 1217 1492
## y1[5] 19.63 19.58 4.43 4.06 12.46 26.72 1.00 1611 1795
## y1[6] 22.68 22.62 4.60 4.14 15.49 30.21 1.00 1924 1933
## y1[7] 25.71 25.72 4.53 4.09 18.42 33.08 1.00 1792 1752
## y1[8] 28.79 28.64 4.67 4.43 21.54 36.39 1.00 1486 1585
## y1[9] 31.97 31.62 4.82 4.37 24.78 40.29 1.00 1767 1903
## y1[10] 34.96 34.76 5.12 4.51 27.01 43.68 1.00 1065 1042
y0=mcmc$draws('y0')
smy0=summary(y0)
grid.arrange(
qplot(xya$y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))
par(mfrow=c(1,2))
hist(xya$y-smy0$median,xlab='obs.-prd.',main='residual')
density(xya$y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')
m1=mcmc$draws('m1')
smm1=summary(m1)
y1=mcmc$draws('y1')
smy1=summary(y1)
xy=tibble(x=x1,m=smm1$median,ml5=smm1$q5,mu5=smm1$q95,yl5=smy1$q5,yu5=smy1$q95)
qplot(x,y,col=I('red'))+
geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
geom_line(aes(x=x,y=m),xy,col='black')
sensitivity/specificity analysis
ex14.stan
estimating sens and spec
data {
int N;
array[N] int x;
array[N] int y;
}
parameters {
real<lower=0,upper=1> p;
real<lower=0,upper=1> se;
real<lower=0,upper=1> sp;
}
model {
p~uniform(0,1);
se~uniform(0,1);
sp~uniform(0,1);
for (i in 1:N) {
y[i]~bernoulli(x[i]*se+(1-x[i])*(1-sp));
}
}
generated quantities {
real ppv;
real npv;
ppv=se*p/((se*p)+((1-p)*(1-sp)));
npv=(1-p)*sp/(((1-p)*sp)+(p*(1-se)));
}n=20
data=list(N=n,
x=sample(0:1,n,replace=T),
y=sample(0:1,n,replace=T))
mdl=cmdstan_model('./ex14.stan')
mle=mdl$optimize(data=data)## Initial log joint probability = -12.8506
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 6 -11.2128 0.000179292 5.9677e-06 1 1 9
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -11.21
## p 0.45
## se 0.78
## sp 0.73
## ppv 0.70
## npv 0.80mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.seeMCMC(mcmc,ch=F)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -17.54 -17.20 1.33 1.15 -19.99 -16.05 1.00 767 828
## p 0.51 0.51 0.30 0.38 0.05 0.96 1.00 460 495
## se 0.73 0.75 0.13 0.12 0.49 0.92 1.00 1918 1304
## sp 0.69 0.71 0.12 0.13 0.47 0.87 1.00 3189 1493
## ppv 0.65 0.73 0.28 0.30 0.10 0.98 1.00 484 570
## npv 0.65 0.73 0.29 0.31 0.10 0.99 1.00 478 617
bimodal distribution, mixed normal distribution
n=100
y0=rnorm(n,0,1)
y1=rnorm(n,-5,1)
y2=rnorm(n,5,1)
y=sample(c(y0,y1,y2),n)
density(y) |> plot()
EM algorithm
library(mclust)
rst=Mclust(y)
summary(rst)## ----------------------------------------------------
## Gaussian finite mixture model fitted by EM algorithm
## ----------------------------------------------------
##
## Mclust E (univariate, equal variance) model with 3 components:
##
## log-likelihood n df BIC ICL
## -254 100 6 -535 -538
##
## Clustering table:
## 1 2 3
## 32 38 30rst$parameters## $pro
## [1] 0.320 0.384 0.295
##
## $mean
## 1 2 3
## -5.1792 0.0431 4.8765
##
## $variance
## $variance$modelName
## [1] "E"
##
## $variance$d
## [1] 1
##
## $variance$G
## [1] 3
##
## $variance$sigmasq
## [1] 1.14
##
##
## $Vinv
## NULLplot(rst)
ex16-1.stan
data {
int N;
vector[N] y;;
}
parameters {
simplex[3] h; //ratio of mix
real m1;
real m2;
real m3;
real<lower=0> s1;
real<lower=0> s2;
real<lower=0> s3;
}
model {
s1~cauchy(0,5);
s2~cauchy(0,5);
s3~cauchy(0,5);
for (i in 1:N) {
vector[3] p;
p[1]=log(h[1]) + normal_lpdf(y[i] | m1, s1);
p[2]=log(h[2]) + normal_lpdf(y[i] | m2, s2);
p[3]=log(h[3]) + normal_lpdf(y[i] | m3, s3);
target+=log_sum_exp(p);
}
}mdl=cmdstan_model('./ex16-1.stan')
data=list(N=n,y=y)
mle=mdl$optimize(data=data)## Initial log joint probability = -635.698
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 99 -265.411 0.000104163 0.0321676 0.3444 0.03444 130
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 100 -265.411 0.000257792 0.0115348 1 1 131
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -265.41
## h[1] 0.70
## h[2] 0.25
## h[3] 0.05
## m1 -2.18
## m2 5.05
## m3 0.71
## s1 3.10
## s2 0.88
## s3 0.04mcmc=goMCMC(mdl,data,smp=2000)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 2100 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 2100 [ 4%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 2100 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 2100 [ 4%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 2100 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 2100 [ 4%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 2100 [ 0%] (Warmup)
## Chain 1 Iteration: 1100 / 2100 [ 52%] (Sampling)
## Chain 2 Iteration: 1100 / 2100 [ 52%] (Sampling)
## Chain 3 Iteration: 1100 / 2100 [ 52%] (Sampling)
## Chain 1 Iteration: 2100 / 2100 [100%] (Sampling)
## Chain 1 finished in 1.3 seconds.
## Chain 2 Iteration: 2100 / 2100 [100%] (Sampling)
## Chain 3 Iteration: 2100 / 2100 [100%] (Sampling)
## Chain 2 finished in 1.3 seconds.
## Chain 3 finished in 1.3 seconds.
## Chain 4 Iteration: 101 / 2100 [ 4%] (Sampling)
## Chain 4 Iteration: 1100 / 2100 [ 52%] (Sampling)
## Chain 4 Iteration: 2100 / 2100 [100%] (Sampling)
## Chain 4 finished in 2.2 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 1.5 seconds.
## Total execution time: 2.4 seconds.seeMCMC(mcmc,ch=F)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -261.34 -260.89 2.58 2.08 -265.73 -258.30 1.00 2237 1424
## h[1] 0.34 0.34 0.07 0.07 0.24 0.46 1.25 11 94
## h[2] 0.32 0.32 0.07 0.06 0.23 0.44 1.24 11 40
## h[3] 0.33 0.33 0.07 0.06 0.24 0.44 1.17 15 72
## m1 -0.06 0.08 3.71 6.19 -5.34 5.09 2.33 4 26
## m2 1.14 2.49 4.28 3.71 -5.35 5.18 2.10 5 29
## m3 -1.34 -2.31 4.16 4.21 -5.47 5.07 2.10 5 26
## s1 1.17 1.12 0.66 0.22 0.82 1.62 1.05 55 3059
## s2 1.15 1.10 0.40 0.22 0.81 1.61 1.06 42 1540
## s3 1.19 1.15 0.29 0.22 0.84 1.67 1.06 39 195
